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Abstract 

We introduce a new discretization scheme for Anisotropic 
Diffusion, AD-LBR, on two and three dimensional carte- 
sian grids. The main features of this scheme is that it is 
non-negative, and has a stencil size bounded by 6 in 2D, 
by 14 in 3D, despite allowing diffusion tensors of arbitrary 
anisotropy. It also has good spectral properties, which 
permits larger time steps and avoids e.g. chessboard arti- 
facts. 

AD-LBR relies on Lattice Basis Reduction, a tool from 
discrete mathematics which has recently shown its rele- 
vance for the discretization on grids of strongly anisotropic 
Partial Differential Equations [3 [7] . We prove that AD- 
LBR is in 2D asymptotically equivalent to a finite element 
discretization on an anisotropic Delaunay triangulation, a 
procedure more involved and computationally expensive. 
Our scheme thus benefits from the theoretical guarantees 
of this procedure, for a fraction of its cost. Numerical 
experiments in 2D and 3D illustrate our results. 



Keywords: Anisotropic Diffusion • Non-Negative Nu- 
merical Scheme • Lattice Basis Reduction. 

We consider throughout this paper a bounded smooth 
domain J7 C K'^, where d G {2, 3} denotes the dimension, 
equipped with a continuous diffusion tensor D. We do 
not impose any bound on the diffusion tensor anisotropy, 
and we are in fact interested in pronounced, non axis- 
aligned anisotropies. Anisotropic diffusion is here under- 
stood in the sense of p3]: the diffusion tensor D(z), at a 
point z £ n, is a symmetric positive definite matrix which 
eigenvalues may have different orders of magnitude. Our 
results are not relevant for isotropic diffusion with a vari- 
able scalar coefficient, as in the pioneering work of Perona 
and Malik 

We address the discretization of the following energy £, 
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defined for u e H^{n): 



We denote 1 1 el 
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|Vu(z)||d(2)'^^- 



v/(^;Miy, for any e 



(1) 



and any 



M in the set S'^ of symmetric positive definite d x d ma- 
trices. Gradient descent for the energy ([!]) has the form 
of a parabolic PDE: 



dtu = div(D Vu) 



(2) 



This equation. Anisotropic Diffusion, is with its variants 
at the foundation of powerful image processing techniques. 
Some variants include curvature terms [1^, or diffusion- 
reaction terms |^. Time varying and solution dependent 
diffusion tensors can also be considered. A general ex- 
position can be found in |14| . where various choices for 
the definition of the diffusion tensor D from the image 
u, adapted to various applications, are proposed and dis- 
cussed. 

Our contribution in the discretization of the energy ([T]) 
results in improved numerical solutions of ([2]), in terms 
accuracy and stability, for a minor increase in complexity. 
This extends to applications, such as Coherence Enhanc- 
ing Diffusion, see the numerical experiments in S|4j which 
involve solving ^ using a solution dependent diffusion 
tensor D = D^; regardless of the fact that the resulting 
nonlinear PDE dfU = div(D„ Vw) may not anymore be 
the gradient descent of an energy. In these applications, 
the diffusion tensor D„ is typically defined in terms of the 
structure tensor |14) of m, in such way that diffusion is 
pronounced within image homogeneous regions, and tan- 
gentially along image edges, but not across edges. 

The 2-dimensional anisotropic diffusion is discussed in 
[14j (chapter 3), where it is proved that a nonnegative 
scheme exists for a stencil of size {2m + 1) x (2m + 1), 
where m depends on the anisotropy of the diffusion ten- 
sor. On the other hand, our approach provides a stencil 
composed of 6 points with non-negative weights. The dis- 
tance from these points to the stencil center depends on 
the anisotropy of the tensor, see Remark [l] 

Consider a scale parameter /i > 0, and a sampling fJ/j 
of the domain Q. on the cartesian grid Z'', rescaled by h: 
with obvious notations 
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We introduce a novel discretization of the energy ([!]), re- 
ferred to as AD-LBR (Anisotropic Diffusion using Lattice 
Basis Reduction). It is a sum of squared differences of a 
discrete map u e L'^{^h) 



£h{u) := h"-^ J2 



E 

zerih eev{z) 



7^(e)|u(z + /ie)-u(z)|2 (3) 



The stencils V{z) C Z'^, z e fi/i, are symmetric and have 
cardinality at most 6 in 2D, 14 in 3D. The coefhcients 
72(e) > are non-negative. They are constructed using a 
classical tool from discrete mathematics, Lattice Basis Re- 
duction, which allows to cheaply build very efhcient sten- 
cils for grid discretizations of Partial Differential Equa- 
tions (PDEs) involving strongly anisotropic diffusion ten- 
sors or Riemannian metrics. This approach was applied 
to anisotropic static Hamilton Jacobi PDEs in [6l , re- 
sulting in a new numerical scheme: Fast Marching using 
Lattice Basis Reduction (FM-LBR). Substantial improve- 
ments were obtained in comparison with earlier methods, 
in terms of both accuracy and complexity. 

The paper is organized as follows. We describe the sten- 
cils of the two dimensional AD-LBR in SjT] and state our 
main 2D result: the asymptotic equivalence of the AD- 
LBR with a finite element discretization on an Anisotropic 
Delaunay Triangulation. Section f|2]is devoted to the con- 
struction of the three dimensional stencils of the AD-LBR, 
and the proof of the cocfhcicnts non-negativity. The more 
technical S|3] details the proof of the 2D equivalence result 
stated in S|T] Two and three dimensional numerical experi- 
ments are presented in fj4j including qualitative and quan- 
titative comparisons with four other numerical schemes. 

1 Description of the scheme, and 
main results 

Our numerical scheme. Anisotropic Diffusion using Lattice 
Basis Reduction (AD-LBR), involves the construction of 
stencils which geometry is tailored after the local diffu- 
sion tensor. This diffusion tensor D is meant to measure 
gradients, as in ([T]). We shall measure vectors through its 
inverse M, which we regard as a Riemanniar{^ metric: for 
all z e fi 

M(z) = D(z)~\ 

An essential feature of AD-LBR is its non-negativity: 
the discrete energy £h{u) is written as a sum (|3| of 
squared differences of values of u, with non-negative 
weights 72(e) > 0. This discretization is consistent if for 
each z € ^1^, and any smooth u, 

iz) E lAe){Vuiz),her. (4) 

eeV{z) 



h'\\Vu{z)\\l 



^The Laplace Beltrami operator associated to M docs not coin- 
cide with div(D V-), unless D is identically of determinant 1. This 
is not an issue for our application. 





Figure 1: Right : the stencils associated to three matrices 
M of anisotropy ratios k{M) equal to 1.1, 3.5, 8 respec- 
tively. The ellipses {z G M^; \\z\\m = 1} are shown left; 
their principal axis is aligned with (cos(7r/3), sin(7r/3)). 
More stencils are shown in [5]. 



Indeed, the left hand side approximates the contribution 
of the "voxel" z + [-h/2, h/2]'^ to the integral Q, while 
the right hand side is obtained by injecting the first order 
approximation u{z 4- he) ~ u(z) + (Vu(z), he) in ([s]). The 
identity Q is in turn equivalent to 



e<£V{z) 



ejee 



(5) 



The next lemma shows how to obtain such a decomposi- 
tion in 2D. We denote by := (— &, a) the rotation of a 
vector u — (a, 6) S by tt/2, in such way that for all 
V G E^: 

(u^, w) — det(u, v). 

Lemma 1. Let ei, 62, es G be such that ei + 62 + ea = 
0, and that a :— det(ei,e2) is non-zero. Then for any 
M€S+: 

M =2Y. 7.e.e, , w^th 7. := 2^2det(M) ' 

l<i<3 ^ ' 

Proof. Note that a — det(e2, ea) — det(e3, ei). 

We first assume that M is the identity matrix Id. De- 
noting D := 2J2'i=i li^i^I 1 we obtain 

{ei,Dei)^2 ^ 7,(e„e^)2 = 2 ^ 7,det(e„ei)2 

l<i<3 l<i<3 

= 20:^(72 +73) = -(ei,e2 -l-ea) = ||ei|p. 



Thus {ei,Dei'^ 



=-L||2 



Likewise (e^,Z)e^^ 



=-L||2 



and (e^ + e^, i^(e^ + e^)) = {ei,Dei) = ||e^f = ||e^ + 
e^lp. Since (ej'-,e^) is a basis of M^, it follows that D = 
Id, which establishes 

In the case M ^ Id, we obtain ^ by applying the above 
to e- := Msei, and M' := Id. □ 

The AD-LBR is based on decompositions ([s]) similar to 
the previous lemma, but with non-negative weights. This 
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construction is discussed below in the 2D case, see ^|2]for 
the 3D case. For that purpose we introduce a classical 
tool of discrete geometry: lattice basis reduction, see [5] 
and references therein. A basis of the lattice 1? is a pair 
(e, /) of elements of 1? , which satisfies 

|det(e,/)| = 1. (7) 

If (e, /) is a basis of Z^, then any g £ 1? can be expressed 
as a linear combination g = ae-\- f3f, with integer coeffi- 
cients a, P G Z. Note also that the coordinates of e are 
co-prime: if e = (61,62) then gcd(6i,e2) = 1, and like- 
wise for /. Our discretization scheme involves privileged 
reduced bases [S], associated to the local metric. 

Consider a fixed symmetric positive definite matrix AI e 
82- The first Minkowski minimum Ai(Af) is the smallest 
norm || • ||m of a non-zero vector in the grid Z^: 

Ai(A/) := min ||e||A/. (8) 

eeZ2\{0} 

Denoting by e a minimizer in the above expression, the 
second Minkowski minimum is the smallest norm 1 1 • 1 1 a/ of 
an element of non-coUinear to e: 

A2(Af):= min |1/|1m. (9) 

A M-reduced basis is a pair (e, /) of minimizers of ([s]) 
and Q . Simple arguments of linear algebra [S] show that 
the value A2(M) does not depend on the choice of the 
minimizer e of ([s]), and that a M-reduced basis is auto- 
matically a basis of Z^, in the sense of Q. We emphasize 
that obtaining a Af -reduced basis, i.e. solving the mini- 
mization problems ^ and (|9|, is both simple and cheap 
numerically. This is the object of Gauss's algorithm [5]: 
initialize (6,/) as the canonical basis of Z^, and 

Do (e, /) := (/, 6-Round((e,Af/)/||/||2,)/), (10) 

The number of iterations is ©(In k(A'/)), logarithmic in 
the anisotropy ratio k{M) of the matrix A/, which is de- 
fined by 

^(M) := V\m\\\M-m = ,max . 

Il"ll=ll'"ll=i II^IIm 

The elements e, /, of a Af-reduced basis, are heuristi- 
cally never very far from being orthogonal, with respect 
to the scalar product (•,A/-). Indeed, since / minimizes 
II • \\m among elements of Z^ non collinear to e, we have 
II/IIm < ||/-6||m and ||/||m < ||/ + e||M- Squaring these 
inequalities, we obtain 

2|(e,A//)|<||e||i,. (11) 



Assume that (e, Mf) < 0, up to replacing / with — /, and 
define g := — e — /. We obtain 

(/,A.fg) = -{e,Mf) - \\f\\l < -^\\f\\l, < 0, (12) 
{g,Me) - -(6, A//) - ||e||2, < -^||e||2, < 0. 

If A/ — M(2;), then we define the weights jz in the 
AD-LBR ^, by 

7.(±6) :=-(/,Afg)/(2detM), (13) 

and likewise permuting circularly the roles of e,f,g. The 
map 72 : Z'^ — ^ IR+ is non-negative, even, identically zero 
except at e, /, perhaps and their opposites. See Figure 
[1] We show in lemma[8]that jz is independent of the choice 
of A/-reduced basis {e,f). In view of Lemma [l] we just 
constructed the decomposition ^ of D(z) = M(z)^^ = 
A/^^, with non-negative coefficients as desired. 

A 3D extension of this construction is proposed in ^ 
The proof is partly computer assisted, as it involves solv- 
ing a 21 X 21 linear system (explicit and with integer en- 
tries) which solution is checked to be non-negative. The 
construction has a similar cost 0{\n k{M)). 

The AD-LBR energy £h : L^i^h) M+, see is 
written in terms of the above constructed coefficients 72, 
and of stencils V{z) C Z^, defined as follows. For each 
z £ Q, we set 

V{z):={eeZ^-jz{e)^0}, (14) 

in the case of periodic or reflected boundary conditions on 
a rectangular domain, and in the case of Dirichlet bound- 
ary conditions on a general domain (extending u by zero 
outside ilh)- In the case of Neumann boundary conditions 
on a general domain, one should set 

V{z; h) := {e G Z'^; 72(e) and z + he e ^h}. 

We have so far established three strongpoints of the AD- 
LBR: 

Non- negativity. Off diagonal coefficients of the symmet- 
ric semi-definite N x N matrix, N — #(f2/i), associ- 
ated to the energy £h sjie non-positive, while diagonal 
coefficients are positive. 

Sparsity. Stencil size is uniformly bounded, indepen- 
dently of the anisotropy of the diffusion tensor: 
#(y(z)) < 6 in 2D, and #{V{z)) < 14 in 3D, for 
all 2 e O. 

Complexity. The construction of the stencil V{z), and of 
the associated coefficients 72, has a logarithmic cost 
©(In k(D(z))) in the anisotropy ratio of the diffusion 
tensor. 
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The next result, Theorem [T] restricted to the two di- 
mensional case, establishes that AD-LBR is asymptot- 
ically equivalent to a more involved and computation- 
ally intensive procedure: a finite element discretization 
of the energy Q, on an Anisotropic Delaunay Triangula- 
tion (ADT, see [Sj and below) of the domain Vl. Under 
the assumptions of this theorem, AD-LBR benefits from 
two additional guarantees, that we state informally and 
without proof. 

No chessboard artifacts. Some numerical schemes for 
anisotropic diffusion suffer from chessboard artifacts, 
in the sense that periodic artifacts develop at the pixel 
level. Such artifacts cannot develop in finite element 
discretizations, since they would lead to high fre- 
quency oscillations of the finite element interpolant. 



and therefore to an increase of the energy (15 1. The 




asymptotic equivalence of the AD-LBR with a finite 
element discretization also rules out these defects. 



Figure 2: Construction of an Anisotropic Delaunay Tri- 
angulation. Left; the grid points p G ri/t, and the Voronoi 
cells Vor;i(p). Right: the triangulation Tin obtained 
(generically) by connecting grid points which Voronoi cells 
intersect. The stencil Vh{p) is represented by arrows. (In 
order to handle the non generic cases where four or more 
Voronoi regions intersect at the same point, an intermedi- 
ate polygonization Qh is introduced in the text.) 



Spectral correctness. The n-th smallest eigenvalue 
\n{h) of the symmetric matrix associated to h~'^£h 
([3]), converges as /i — towards the n-th smallest 
eigenvalue A„ of the continuous operator — div(D V), 
for any given integer n > 0. This follows from a simi- 
lar property of the finite element energy £'i^ ( [T5| , and 
from the asymptotic equivalence (|16|). 



Our convergence result. Theorem [T] below, is special- 
ized to the case of a square periodic domain, which covers 
reflecting boundary conditions frequently used in image 
processing. Since the grid discretization must be compat- 
ible with the boundary conditions, any scale parameter h 
appearing in the rest of the paper is assumed to be the 
inverse of a positive integer: 

h e {1/n; n > 1}. 

Theorem 1. Let i7 be the unit square [0, , equipped with 
periodic boundary conditions. Let D : 11 — > 5*^ be a (pe- 
riodic) diffusion tensor with Lipschitz regularity, and let 
M := D^^. When h is sufficiently small, the periodic Rie- 
mannian domain (17, M) admits an Anisotropic Delaunay 
Triangulation Th, with collection of vertices il^ := Q,C\hl? . 
For u G i^(il/i), define 



|V(Ir, w)(2;)||d(2)0^z, 



(15) 



where I7- denotes the piecewise linear interpolation op- 
erator on a triangulation T ■ Then for some constant 
c = c(D), independent of u and h. 



(1 - ch)£h{u) < < (1 + ch)£h{u). 



(16) 



The proof of this result is postponed to SjSj but for 
the sake of concreteness, we describe the concept of 



Anisotropic Delaunay Triangulation (ADT) [5j. In the 
rest of this introduction, and in ^|3j we assume as in Theo- 
rem [T] that the diffusion tensor D is defined on the square 
[0, 1]^ and satisfies periodic boundary conditions. We ex- 
tend it, as well as its inverse the metric M, to the whole 
plane by periodicity. 

We specialize the concept of ADT 5J , to the domain 
and the collection of vertices hZ'^. For that purpose, we 
introduce some notations. For all p, g G K^, we denote by 
Sp{q) the distance from p to q, as measured by the metric 
at the point p: 



Spii) ■= \\q-p\ 



M(p)- 



We denote by Ah{q) the least distance from a point q G M^, 
to the grid hZ'^: 



Ah{q) ■■= min 5p{q). 



(17) 



We introduce the Voronoi cell Vor/i(p) of a grid point p G 
hZ"^ , which is the collection of points q G closer to p 
than to any other grid point: 

Movhip) ■■= {q e 5p{q) = A„(g)}. (18) 

The collection of Voronoi cells is referred to as the Voronoi 
diagram, see Figure [2] (left) . A Voronoi vertex is a point 

q e 

intersect: {yoTh{Pi))^^i, k > 3, p. 
a dual Voronoi cell Tg, defined as the convex hull of the 
points (ft)*Li- 

The geometric dual Qh, of the Voronoi diagram, is de- 
fined as the collection of all dual Voronoi cells Tg. Note 
that, generically on the metric M, no more than three 
Voronoi regions can intersect at any point in M^, thus 



^ at which at least three distinct Voronoi regions 

G hZ^. We attach to q 
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the elements of Qh are generically triangles. If h is small 
enough, we show in Sjs] (using the Dual Triangulation The- 
orem in that Tg is a strictly convex polygon, of vertices 
iPi)i=i with the above notations, and that Qh is a polygo- 
nization (generically a triangulation) of , with vertices 

Since the metric M and the vertices h'E'^ are periodic 
(recall that h — 1/n for some integer n > 1), arbitrarily 
triangulating the elements of Q/i, respecting periodicity, 
yields a periodic triangulation Th- 

Definition 1 (ADT, Labelle and Shewchuk 0). The tri- 
angulation Th obtained by the above construction is re- 
ferred to as an ADT of the domain M^, with collection 
of vertices hlP' , and underlying Riemannian metric M. 
Since Th is iP' -periodic, we also regard it as an ADT of 
the periodic unit square Q.. 

We establish in §3.1 the existence of the ADT Th- Sub- 
section §3.2 is devoted to the study of reduced bases: 
their characterization, uniqueness and stability properties. 
We study in §3.3 the finite element stencils, defined for 
p e hi? by 

Vh{p) ■■= {e e I?] [p,p + he] is an edge of %}■ (19) 

See Figure [2] (right). We show that Vh{p) coincides with 
the AD-LBR stencil V{p), unless the lattice 1? admits a 
basis almost orthogonal with respect to the scalar product 
associated to M(p), see Lemma 10 This is tied to the 



fact that orthogonal grids admit several (usual) Delaunay 
triangulations. Overcoming this technical difficulty, we 
conclude the proof of Theorem [ij 

Note that the construction of the ADT Th is not easy 
to parallelize, in particular when anisotropy is pronounced 
since the Voronoi' regions of far away points interact. The 
construction of Th also involves solving polynomial equa- 
tions of degree four, because Voronoi regions boundaries 
are conies, and Voronoi vertices must be identified at their 
intersections. In contrast, the stencils of the AD-LBR are 
independent of each other, and their construction is essen- 



tially contained in the two lines program ( 10 ) 



2 Three dimensional stencils 

We construct in this section the 3D stencils of the AD- 
LBR. We also discuss in Remark [T] the euclidean radius of 
these stencils, in 2D and 3D. 

We denote by eiZ + • • • + e^Z the sub-lattice of Z'' gen- 
erated by vectors ei, • • • , e^. G Z''. This sub-lattice equals 
{0} by convention if fc = 0. We say that (ei, • • • , e^) G 
{U^Y is a basis of Z'' if eiZ + ■ • • + e^jZ = l/, or equiva- 
lently if 

I det(ei, • • • ,ed)\ = 1. 






Figure 3: Right: stencil of the AD-LBR, for a symmetric 
matrix of eigenvector M of anisotropy ratio k{M) equal to 
2 (top) or 6 (bottom). The anisotropy is of "needle" type: 
the two largest eigenvalues of M are equal, and the needle 
orientation is given by the vector (4,2,3). The ellipsoid 
{z e K^. ii^ii < X} is shown left. 



Definition 2. Given M E S'^ , d < 4, we say that a basis 
(ei, • • • , Cd) ofU^ is a M-reduced basis if for all 1 < i < d 

\\ei\\M e argmin{||e||M; e G Z"* \ (eiZ H h e^^iZ)}. 

There always exists a M-reduced basis in dimension d < 
4, see [5]. The norms 



A,(M) := ||e,||M, 



(20) 



of the elements {ei)f^i of a M-reduced basis, are called the 
Minkowski minima, and are independent of the choice of 
M-reduced basis. Such bases can be obtained via a gen- 



eralization of Gauss's algorithm (10), still of complexity 



0{\nK{M)), see [9] for the proof, and Theorem 1.5 in [6] 
for a discussion. In dimension d > 5 there may not exist 
any M-reduced basis in the sense of Definition [2] and the 
relevant notion is Minkowski reduction [S]. 

The AD-LBR, three dimensional, 14 element symmetric 
stencil V{z) is obtained by symmetrizing the 7 element 
stencil V described in the next theorem, for M := M(z). 
This stencil was originally designed for solving 3D static 
Hamilton Jacobi (HJ) PDEs in [6,. It is not yet clear 
wether the 4D stencil proposed in [J for HJ PDEs, can be 
used to extend the AD-LBR to four dimensions. 



Theorem 2. Let M £ 

M-reduced basis. Let u 



, and let (61,62,63) be a 
:= ±6(j(2); and 



±e, 
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w = ±eo-(3), where the permutation a of {1,2,3} is such 
that \ {v,Mw)\ > \ {u,Mv)\ > \ {u,Mw)\, and the signs are 
chosen so that {u, Mv) > and {u, Mw) > 0. Thus 



— 1. We denoted by (ai)i=i ^-^d (/3i)f=i the coefficients of 
a and /3. The range of admissible values of t is thus 

max{0, max{— ttj-; /3j = —1}} l£ t < min{ai; /3i — 1}. 



|(t;,Mu;)| > {u,Mv) > {u,Mw) > 0. 



(21) 



If this range is non-empty, then 



> whenever 



Then there exists a 7 element family V = (ei)i=i; <^i^d 
non-negative coefficients (tOI^i; such that 



CiCA 



l<i<7 



and 



(n) If {v,Mw) < 0, then 

V := (u, V, w, u — V, u — w, V + u), u — V — w). 

(p) If{v,Mw) > 0, then 

V :— {u, V, w, u ~ V, u — w, V — w, u — V + w). 

Proof. Let A be the 3x3 matrix which columns are 
{u,v,w). We assume, up to replacing M with A^MA, 
that {u, V, w) is the canonical basis of K'^. 

We denote hy u ■ v the scalar product {u, Mv), and by 
the squared norm = {u, Mu) = u-u. Recall that 

the inverse of M is positively proportional to its coma- 
trix M := det(M)M-^ The coefficients {M,j)lj^-^ of the 

comatrix are 2 x 2 minors of M, hence quadratic forms over the six variables {u^ 



/3i(3j — —1, and > whenever /3j — 1. Conversely, if 
these inequalities are satisfied, then the range of admissi- 
ble t is non-empty, and Theorem [2] is proved. 

We thus only have to prove the non-negativity of a finite 
number, twelve in each case to be exact, of quadratic forms 
{ai, or ai + aj, for the values of i,j given above) in the 
quantities (u^, v'^ , w'^ , u ■ v, u ■ w, v ■ w). We know that 
the six following linear forms of these quantities are non- 
negative: 

• In case (n): u ■ w, u ■ v — u ■ w, —v ■ w — u ■ v. 



2u ■ V, v + 2v ■ w, w +2v ■ w. 



In case (p); u ■ w, u ■ v — u ■ w, v ■ w — u ■ v. 



and 



2u ■ V, v — 2v ■ w, 



2v ■ w. 



In each case, the three first follow from ( 21 ), and the three 



last are obtained as in ( 11 ) (use the fact that ||ei±ej ||m > 
IjeillM for any distinct Multiplying these six lin- 

ear forms by one another, we obtain 21 distinct quadratic 
forms which take non-negative values, and happen to be a 
basis of the 21-dimensional vector space of quadratic forms 



in its coefficients (Mij)?^-^j^. For instance, recalling that 
{u, V, w) is here assumed to be the canonical basis of M.^: 

Mil = M22M33 - M23M32 = v^w^ - {v ■ w)^. 

M21 = M32M13 - M12M33 ^ {v ■ w){u ■ w) ~ {u ■ v)w'^. 

The rank 1 symmetric matrices zz^ , where z ranges over 
the six first elements of V, are easily checked to be a basis 
of six dimensional space of 3 x 3 symmetric matrices. The 
vector a € of the coefficients of M on this basis is given 
below in case (n): 

a := (Mil + M12 + Mi3, Af2i + M22 - M23, 

M31 - A#32 + M33, -M12, -Mi3, A#23)- 

In case (p), replace M23 with — M23 (and Af32 with 
— M32) in the above expression of a. The rank one 
matrix zz^ , where z is the last element of the stencil 
V, can also be decomposed on this basis, with coeffi- 
cient vector j3 := (— 1, — 1, — 1, 1, 1, 1) in case (n) (resp. 
/3 := (1,-1, 1,1, -1,1) incase (p)). 

As a result, we can write the comatrix M as the linear 
sum of zz"^ , where z ranges over V, with coefficient vector 
{a-tj3, t) e M"^, for any t £ R. The proof of Theorem[2]will 
be complete when we find t e , such that a~t/3E M.^ . 

Such a value of t will satisfy t < ai for each 1 < i < 6 
such that (3i = 1, and t > —aj for each j such that (3j = 



V , w , u ■ v,u ■ w,v ■ w). It also 
happens that the coefficients of the 12 of quadratic forms 
interest (a^ or ai + aj , see above) on this basis are non- 
negative; this is shown by inverting a 21 x 21 explicit ma- 
trix with integer entries, see the ancillary file. Hence they 
take non-negative values, which concludes the proof. □ 

Remark 1 (Stencil radius). Consider a matrix M G S'^ , 
where d G {2, 3}, and the associated stencil V used in 
the AD-LBR. As discussed in the introduction, the stencil 
cardinality ^{V) has a direct impact on the scheme com- 
plexity, and is hopefully bounded independently of M: by 
6 in 2D, and by I4 in 3D. The stencil euclidean radius 

r := max{||w|| ; v G V} 

can however also have an indirect impact on computation 
time. A large radius may e.g. provoke non-local memory 
accesses and cache misses, which can raise crucial issues 
in the prospect of a GPU implementation. 
Observe that, for any e G M.'^ 



\M- 



2 < 



< 



|M|| 



Thus, in view of the stencil construction, 
r < d\d{M)\\M-^\\^ . 



(22) 



(23) 



Inequality ( 22 ) also implies 



\M-^\\-^ < Xi{M) < ■■■ < Xd{M) < ||Af||5, (24) 
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see e.g. Proposition 1.6 in fB^. Combining ( |23[ ) and (24), 
we obtain the worst case upper bound r < dn^M). 

An average case estimate of is proved in Q/. Denot- 
ing by Od the group of rotations o/M'', equipped with the 
Haar measure: 



\d{R^MR)dR < Cddet{M)-- 



(25) 



Let r be the average value, over grid orientations, of the 
stencil radius. Combining (23) and (|25[) w e obtain, in two 
dimensions, the estimate r < 26*2 \/ k{M). 

3 Equivalence to a finite element 
discretization 

This section is devoted to the proof Theorem [l} the 
asymptotic equivalence of AD-LBR with a finite element 
discretization. We use the notations of fjl] The existence 
of the ADT Th is established in the first subsection, for h 
sufficiently small, as well as a few of its properties. The 
second subsection is devoted to the study of M-reduced 
bases. Theorem [l] is proved in the third subsection, by 
comparing the stencils of the AD-LBR and of the finite 
element discretization. 



3.1 Existence of an ADT 

We introduce three positive reals i>i, 1/2,1^^, defined by 

:= min{|| M(z)"^||" 
max{|| M(z 



2 • z e 



(26) 



V2_ 



Thus for all p, r e 

vAv - ''ll < ^p{r) := \\p - r||M(p) < 1^211^ - ^11- (27) 
Lemma 2. • For all r £ M^, one has A/j(r) < i'2h. 

• If p,q £ hi?, and r £ YoTh{p) H Vor/i(q), then 
\\p — ^11 < v^h and ||p — q\\ < 2v^h. 

Proof. First point. Rounding the coordinates of r to a 
nearest multiple of h, we obtain a point p £ hi? such that 
\\p — A\ — ^- Recalling (27) we obtain dp{r) < V2h, and 



therefore A;i(r) < 1/2/1 in view of (17) 



Second point. We have i/i||p — r\\ < dp{r) A;i(r) < 
1^2^. Thus \\p — r\\ < v^h, and likewise \\q — r\\ < v^h. 
Finally, by the triangle inequality, fzH < ||p — ?'|| + \\q — 
r\\ < 2v^h. □ 

Following the notations of [5|, we denote by T{p,q), 
p, g £ M^, the smallest constant r > 1 such that 

r~Mp(r) < dg{r) < T6p{r), for all r £ M^. 



Equivalently, in the sense of symmetric matrices, 

M{p) < M{q) < M{p). (28) 

We also define a quantity > 1, closely related to the 
modulus of continuity of the metric M: 



Th := max{T(p, g); \\p - q\\ < 2v^h}. 



(29) 



One has r/j — ^ 1 as /i — 0, for any continuous metric M 
(indeed M is periodic and therefore uniformly continuous) . 
If M is Lipschitz, as assumed in Theorem [T] then t/j = 
l + 0{h). 

We assume, in the rest of this subsection, that h is suf- 
ficiently small so that 



(30) 



Lemma 3. • If p,q £ hi?, p ^ q, and r £ Voihip) H 
YoThiq), then 6p{r) < 5p{q) / T{p,qY - 1. 

• The geometric dual Qh of the Voronoi diagram is, as 
announced in f|7| a polygonization of into strictly 
convex polygons, with vertices hi? . 

Proof. First point. We may assume that T{p,q) > 1, oth- 
erwise there is nothing to prove. The second point of 
Lemma [2] implies that \\p — q\\ < 2v^h, thus 

^r{p,qY-l < < 

On the other hand 6p{q) > i^i\\q — p\\ > I'lh, and Sp{r) < 
A/i(r) < 1^2^. The announced inequality follows. 

Second point. We apply Theorem 7 (Dual Triangulation 
Theorem) in j6] . Since the domain has no boundary, it 
suffices to check that all the Voronoi arcs and vertices are 
wedged, see [5| . This condition means that for any p,q £ 
hi? such that p ^ q, and any r £ Vor;j(p) n Vor^(g), one 
has (r — q) M(g)(p — g) > 0, and likewise exchanging the 
roles of p and q. Heuristically, it expresses the acuteness 
of some angles measured in the local metric. Lemma 5 in 
[5 shows that this condition follows from the first point 
of this lemma, which concludes the proof. □ 

We recall that Th is the triangulation obtained by arbi- 
trarily triangulating the polygonization Qh of the previous 
lemma, respecting periodicity, see Definition [T] Generi- 
cally Qh is already a triangulation, hence Th — Qh: see 

§1- _ 

We establish in the next lemma a few properties of the 
elements oiTh. Note that the vertices p, g, r of any triangle 
T £ Th satisfy by construction 



Vor,,(p) n Vor,,(g) n Vor,,(r) ^ 0. 



(31) 



The Voronoi regions Vor^i, and the triangulation Th are 
illustrated on Figure [2] 
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Lemma 4. Denote by he, hf, hg the edges of a triangle 
T ^ Th, where e, f,g G 1? are oriented so that e+f+g — 0. 
Then 

Diameter: max{||e||, ||/||, ||5||} < 2;^^. 

Area: |det(e,/)| = 1, thus \T\ = h'^/2. 

Acuteness: (e, M(2;)/) < 6f^, for any z G T, where 9h '■— 
^l(3 + 9T|J(r|,-l). 

Proof. First point. Wc denote by p, q, r the vertices of T, 
ordered in such way that p-\-he = q, q+hf — r, r-\-hg = p. 



The announced estimate follows from (31), and from the 
second point of Lemma [2] 

Second point. Since 7/i is a conforming triangulation, 
the intersection of T with the collection hi? of all vertices 
of Th consists of only three points: the vertices p, q, r of 
T. Thus the triangle of vertices — e, 0, /, homothetic to T, 
contains no point of integer coordinates but its vertices. 
This imphes that (e, /) is a basis of Z^, hence | det(e, /)| = 
1, as announced. 

Third point. The pairwise distances between p,q,r are 
bounded by 2i'^h (first point), and since z e T so are the 
pairwise distances between p, q, r, z. Defining s := p — q + 
r G hi?, and observing that \\s — p\\ — \\r ~ q\\ < 2v^h, 
we find that the pairwise distances between p, q, r, z, s are 
bounded by 4i/_^/i. 

Let X G Voihip) n Vor/j(gf) n Vor;i(r). We have Sp{x) — 
Sq{x) — Sr{x) — Ah{x) < Ss{x), thus 

5,{xf > 5p{x)'^ - 6g{x)^ + 6rixf. (32) 

(For intuition: in a classical Delaunay triangulation, x 



would be the circumcenter of T, and (32) state that s 
is outside the circumcircle.) Denoting M :— M(z), and 
S := Ah{x), we obtain 

\6p{xf^\\x^p\\U<6p{xf(r{p,zf-l) 

<SHtI,,-1), (33) 

using Lemma |2] and likewise for q,r. We also have 

5s{x) = \\p - q + r - x\\m{s) 

< \\P - a:||M(s) + Ik - x\\Mis) + \\r - x\\m{s) 



Thus, proceeding as in (33 1, 
Mx) - \\s - x\\U < Ss{x)'{ri,, - 1) < 9<5V|,,(t|, - 1). 



Injecting in (32) these estimates of Sif{x), * G {p,q,r, s}, 
and using the fact that S < hh'2, see Lemma [2] we obtain 
after expansion the announced estimate of (e, M f) . □ 



We next rewrite the finite element energy ( 15 ) in 
a form that can be easily compared with the AD-LBR 



energy £h ([3]), which is written in terms of coefficients 7p 
and stencils V{p), p G Qh. 

We denote by ipp : M the piecewise linear func- 

tion on Th such that ^Pp{p) — 1, and fpiq) = for any 
vertex q G hZ'^ distinct from p. This is the classical "hat 
function" encountered in finite element analysis. For all 
p G hl'^, e G Z^, we define 

lpie)--=-lf (V(^;j(z),D(z)V(p,V,(z))dz (34) 



Clearly, 7p (e) = if + he] is not an edge of Th, in 
other words if e does not belong to the stencil Vh{p), de- 
fined in (19 1. We express in the next lemma the finite 



element energy (15), in terms of the stencils Vh and of 



the (potentially negative) weights 7^. 

Lemma 5. For any u G L'^{Q,h), extended by periodicity 
to hl"^ , one has 

£ki^)-Y. E ^p(^)l"(P + ^^)-"(P)l'- (35) 

Proof. For any triangle T G Th, and any p,q G hl'^, we 
denote 

ST{p,q) / {\7ip';iz),-D{z)V^'^^iz))dz. 
Jt 

Clearly ST{p,q) = if (7 or p is not a vertex of T. The 
coefficient 7p(e), e G Z^, is thus given by the following 
sum with at most two non-zero terms: 



Tp(e) = -\Y1 ST{p,P + he). 



(36) 



Ten 



Let p,q,r G /iZ^ be the vertices of a triangle T G Th- 
Since the sum (fp+(fg+ is constant on T, equal to 1, 
it has a null gradient on T, and therefore 

st{p,p) + ST{p,q) + ST{p,r) = 0. 

Using this relation, and the two similar ones obtained by 
permuting circularly p, q, r, we obtain 



|V(Ir.«)(^)llD(.)c'^ 

= u{pYst{p,p) + u{qfsT(q,q) + u{r)'^ST{r,r) 
+ 2u{p)u{q)sT{p, q) + 2u{q)u{r)sT{q, r) 
+ 2u{r)u{p)sT{r,p), 

= - st{p, q){u{p) - u{q)f - sriq, r){u{q) - u{r)f 
- ST{r,p){u{r) - u{p)f. 

Summing this expression over all T G Th, and combining it 



with (36), we obtain (35 1, which concludes the proof. □ 
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Finally, we provide an approximation of the coefh- 
cients 7^" which will be easily compared with the AD-LBR 



weights 7p (13). 

Lemma 6. Consider an edge + he] ofTh, shared by 
the two distinct triangles T,T' G Tu- Let hf,hg (resp. 
hf ,hg' ) he the two other vector edges of T (resp. T'), 
oriented so that e + f + g = (resp. e + f + g' = 0). 
Then, with M := M{p) 

where eh '■— 2t^l max{|| D(a;) — D(?/)||; — ?;|| < 2v^h}. 

Proof. We assume, up to exchanging / and g, that [p,p — 
hg] is an edge of T. We denote a := det(e, /) G {—1, 1}, 
see Lemma |4] (Area), and observe that a = det (/,(?) = 
deg{g, e). Let 7 be the constant value of Vip^ on T. Then 
i'jjhe) = —1 and {'j,hg) = 1. These two independent 
linear identities are satisfied by af-^ /h, which is thus equal 
to 7. We have shown that Vtpp = af^ /h on T. 

Denoting q := p + hg, we obtain likewise ViyjJ — ag^ /h 
on T. Let D D(p) = A/~^, and let R be the rotation 
by 7r/2. Then 



2 

]^{f,R'DRg)^\ 



ag 



M 



f5 



{f,Mg) 



detAr/ 2detM 
Therefore, using Lemma |4] (Diameter) in the last step, 

{f,Mg) 



(V^;\D(^)V^,^)dz. 



2det M 



^(V^^, (D(z) - D(p))V<^J)dz 

h'^ 2V^ 2V^ rn ^/ N M, 

- y^-^inax{||D(z)-D(p)|| 



Proceeding likewise on T', and recalling (34) (or (36l), we 
conclude the proof. □ 

3.2 Some properties of M-reduced bases 

We establish some technical properties of M-reduced 
bases, thanks to which we will be able to compare in §3.3| 
the "geometric" construction of the ADT finite element 
stencils V/i, with "algebraic" construction of the AD-LBR 
stencils V . 

Lemma 7. Let M e S^, let ei, • • • ,e„ e 71? , n> \, and 

let £ e { — 1, 1}. Assume that for all 1 < i < n, with the 
convention e„_|_i := ei; 



det(ei,ei+i) = e, 
(e,,Mei+i) > 



1 



min{lle»|lM, |le,+i|lM} . 



(37) 
(38) 



Then any M-reduced basis (e, /) satisfies 



{ej} C {ei. 



i}. Our objective is to show 
that z cannot be the element of a M-reduced basis, and 
we may therefore assume that z has co-prime coordinates. 



Proof Let z € Z2\{ei, ■ 



It follows from (37 1 that the closed polygonal line of 
consecutive vertices ei, • • • , e„, circles at least once around 
the origin. Hence z — aci -\- /3ei+i, for some 1 < i < n 
and some a, /3 > 0. Since | det(ei, ei+i)| = 1, a and (3 are 
integers. Since z ^ {e-i,ei+i\, a. + (3 > 2. Since z has 
co-prime coordinates, a/3 7^ 0. 

Assuming without loss of generality that ||ei||M > 
we obtain using (38) 



2 

M 



a 



2 

M 



-illlf + 2a/3(e, 
a/? min{||ei| 
ll + {a^+P^-^-ame^+l\\lI■ 



> a 

> lie 



Me 



2 

M 



-^^l|ei+i||L 



Ml 



1+1/ 



Observing that a^ + (3^ — I — a/3 > for all a,/3 € [1, oo[, 
we obtain ||z||m > ||e,;||A/- Since and e^+i are linearly 
independent, we have ||ei||A/ > A2(M). Finally ||2;||m > 
A2(M), hence z cannot be the element of a M-reduced 
basis, which concludes the proof. □ 

We use in the following the shorthand 

{ai, • • • ,a„, and opposites} 
:= {air ■ ■ , U {-ai, • • • , -a„}. 

Given M S 5^, and a M-reduced basis (e, /) of Z^, we 
denote fJ.{M) := |(e, M/)|. This value can be expressed 
in terms of the Minkowski minima and thus does not de- 
pend on the particular choice of M-reduced basis. Indeed, 
recalling the identity 

(e,M/)2+det(M)det(e,/)2 = |le|li,,||/||2„ 

we obtain 

^(M) = I (e, Mf) I = v/Ai(M)2A2(M)2-det(M). (39) 

A vanishing value, Ijl{M) — 0, indicates that the lattice 1? 
admits an il/-orthogonal basis. 

We next show that the stencils of the AD-LBR do not 
depend on the choices of reduced bases, as was announced 
in the introduction. 



Lemma 8. The weights 73:2^^ M+, defined in (13) 
and used in the AD-LBR, do not depend on the choice of 
Mlz) -reduced basis. 

Proof. We denote M :— 'M.{z), and consider two M- 
reduced bases (e,/), (e',/'). We assume as in the in- 
troduction that {e,Mf) < and {e',Mf) < 0, up to 
changing / or /' into its opposite. Let g := —e — f and 
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g' := — e' — /'. Our objective is to show that the weights 
7,7' : 1? M+, defined by (13) and associated to the 



bases (e,/), (e',/'), are identical. 

RecaU that {f,Mg), {g,Me), {f,Mg'), {g',Me') are 
negative, see (12). Applying Lemma [t] to the family 



/ f I pl II rl\ 

(e ,~g J ,g ) 

we obtain that 

{e, /} C {e', /', g' , and opposites}. 



(40) 



We used ||/||m < K(Af)||e|| a/, see ([24]), in the fourth line. 
Replacing a and r with their assumed upper bounds, we 
obtain 2(e, M'/) < ||e||^^,. Proceeding likewise, we ob- 
tain 2|(e,M'/)| < min{||e|||,,,||/|l2^,}. We may there- 
fore apply Lemma[7]to M' and (e,/, — e,— /), and obtain 
{e'j /'} C {e, /, — e, — /} as announced. □ 



3.3 Comparison of the stencils 

We assume in this subsection that the scale parameter h 
is sufficiently small. Our assumption is stronger than the 



If fi{M) ^ 0, then (e,M/) and {e',Mf') are negative, 
and not merely non-positive, thus {e, /} C {e',f',g'}, or 
{e,/} C {-e',-f,-g'}. Since e + / + g = = e' + 
f+g', it follows that {e, /, g} = {e' , /', g'}, or {e, /, g} = 
{— e', — /', — 5'}. The expression (13) thus implies 7 = 7', 
as announced. 

If At(M) = 0, then (e, Af/) = = {e',Mf'). Using 
again ([T3| we find 

7(±g) = -(e,M/)/(2detAf)=0 (41) 
7(±e) = ||/||2,/(2detAf), 
7(±/) = ||e||l,/(2detAf), 

and likewise for 7', e', /',<?'. Note also that H^'lllf = 



one used in [3.1 see (30 1, hence in particular there exists 



an Anisotropic Delaunay Triangulation Th- More precisely 
we assume that 



1 



and Oh < 0t 



■— 



(43) 



See ( p6| , ([29|, and Lemma|4]for the definition (i^i, v^), 
and Oh respectively. For Lipschitz metrics, t/j = 1 + 0{h) 
and Oh ^0{h). 

Our objective is to compare the stencils V{p), Vh{p), 
of the AD-LBR ([M]) and of the ADT finite element dis- 
cretization (19) respectively, at a point p € hi?. The 



next lemma shows that they are equal unless the lattice 
1? is almost orthogonal with respect to the local metric; 



!m + II/'IIm > >^2{M)\ hence e and / are different from a property quantified via ^(M(p)), see (l39) 



g' and -5'. It follows from (|40]) that {e, /} = {eie',£2/'} 
for some ei,e2 G {~li 1}- This implies 7 = 7' in view of 
(41), and concludes the proof. □ 



We saw in the introduction ( 11 ) that 
2^l{M) < Ai(Af)^ 



(42) 



The next lemma establishes weak unicity and stability 
properties for Af-reduced bases, when the inequality (42) 
is not saturated. 

Lemma 9. Consider M,M' £ S2 , o, AI -reduced basis 
(e, /), and a M' -reduced basis (e',/'). Let t > I be such 
that T~'^M < M' < t'^M , in the sense of symmetric ma- 
trices. Assume either: 

• 2fi{M) < Ai(Af)2, andr^l (i.e. M' ^ M). 

• 4^(A-f) < \i{Mf, and < 1 + \k{M)-^ . 
Then {e', /'} C {e, /, and opposites} . 

Proof. Denoting a := 2/i(Af)/Ai(Af)^, we obtain: 

A{e,M'f) = \\e + f\\lt.-\\e-f\\lj, 

2 

M 



Lemma 10. Letp e hi?, and let M := M(p). If fi{M) > 
Oh, then Vh{p) = V{p). In any case, one has for any M- 
reduced basis (e, /).' 

Vh{p) D {e, /, and opposites} (44) 
Vh{p) C {e, /, e -f /, e - /, and opposites}. (45) 

Proof. We assume that {e,Mf) < 0, up to replacing / 
with — /. Let T £ Th he a triangle containing p, and let 
hex, he2, he^ be the edges of T, oriented so that ei -\- e2 + 
63 = 0. Using Lemma [4] ( Acuteness) , we obtain for all 
1 < i < 3, with the convention 64 := ei 



{e^,Me 



i+l. 



< 



h<0o< \vl 



< 



Denote E :— {61,62,63}, and —E :— {—61, —62, —63}. 
Applying Lemma [t] to M and the points (61,-63,62, 
—61,63,-62), we obtain that {e,f} C E Li [—E). Up 
to exchanging E with —E, we thus have {e, /} C i? or 
{6, — /} C E. Since the elements of E sum to zero, we 
conclude that 



<A\e + f\\\,-r-'\\e-fr 



E^{e,f, 



'-)i\\e\\lj + \\frM) + 2ir^+r-'){e,Mf) 



< ((r^ - r-2)(l + «(A/)2) + air^ + r-^m4li 
<((r4-l)(l + ^(Af)2)+a(r4 + l))||6fM,. 



/} or i? = {6,-/,-6 + /}, (46) 



which implies ( |45[ ). 

If fJ.{M) = 1(6, Af/)| > 9h, then Lemma [I] (Acuteness) 
forbids the second case in (46). Thus E = {e, /, —6 — /}, 
and therefore Vh{p) C V{p). 
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Figure 4: Consider a point p e hi? ^ and denote M :— 
M(p). From left to right: ellipse {||z||m < 1}, AD-LBR 
stencil V{p), stencils W{p) C Vh{p) C W'{p). For W(p) 
and W'{p) we assumed that ^(M) < (?oi otherwise they 
are equal to V{p). Note that V{p), W{p), W'(p) only 
depend on M, while Vh{p) depends on the structure of 
the triangulation Th- 



Let T G 7^ be a triangle containing p and intersecting 
the half line L := {p + re; r > 0}. We know (46 ) that e is 



a vector edge of T. The corresponding edge segment must 
be [p, p + he] , since otherwise TnL would be empty. Thus 
e G Vh{p). Applying the same argument to — e, /, — /, we 
obtain (|44|. 

If IJ.(M) > 9fi, then e + / is also a vector edge of any 
triangle T d Th containing p, since we eliminated the 
second case in (46 1. Reasoning as above we find that 
{e + /, -e - /} C Vh(p), and therefore V{p) C Vh{p). 
Thus V{p) = Vfi{p)- This concludes the proof. □ 

We introduce new stencils W{p), W'{p), where p € M^, 
M := M{p), defined as follows. If ^(M) < Oq, then denot- 
ing by (e, /) a M-reduced basis, 

W{p) :— {e, /, and opposites}, (47) 
W'{p) := {e, f,e + f,e- /, and opposites}. (48) 

On the other hand, if ^{M) > 6^, then 

W{p) := V{p) =: W'{p). (49) 

Note that W{p) C Vh{p) C W'{p), for any p € /iZ^, since 
Oh < Oq by assumption ( 43 ) . 



Lemma 11. The stencils W{p), W'{p), do not depend on 
the choice of M.(p) -reduced basis. 

Proof. Let M := M{p). If fj.{M) > Oq, then W{p), W'{p) 



are defined by (49 1, hence here is nothing to prove. Other- 
wise we obtain /x(M) < 6*0 < J^?/4 < Ai(M)V4. Hence, by 
Lemmajoj any two M-reduced bases (e, /), (e', /'), need to 
satisfy {e', /'} C {e, /, and opposites}. In view of (47) and 
(48), they thus yield the same stencils W{p), W'{p). □ 



W,W': for u € T'^{Q,h), extended to hZ'^ by periodicity, 
^h{u) -.^ ^ ^ \u{z + hg)-u{z)\'^, 

Kiu):^ E E \uiz + hg)-uiz)\'. 
zefifc gew'{z) 

The outline of the proof of theorem [T] is as follows. We 
prove in Lemmas Tsj 14 and 12 respectively that for any 
ueT^iflh): 

14 (u) - £hiu)\ < CoiOh + eh)K{u) (50) 
< CiThiu) (51) 
M^) < C2£h{u), (52) 

where the constants Cq, Ci, C2 only depend on the metric 
M. Combining these inequalities, and recalling that 6h = 
0{h) and Eh = 0{h) for Lipschitz metrics {sh is defined 
in Lemma [6]), we obtain 

\£h{u) ~ £h{'u)\ < ch£h{u), 



for some constant c = c(M). This establishes (16), and 
concludes the proof of Theorem [T] 

For each p G M^, we denote by r?p, 77^ : — )• {0, 1}, the 
characteristic functions of W{p) and W'{p) respectively. 
The proofs of (50) and (52) immediately result from the 



comparison, in Lemmas 13 and 12 respectively, of the co- 
efficients 7p, 7p, 77p, rfp appearing in the expressions of 

£hi£'hi^hT^'h- 

Lemma 12. Tor any p G M^, one has on 1? 

Vp < C2IP, with C2 2df^^ max{det M(g); q e Q}. 
Proof Let M := M(p), and let (e, /) be a M-reduced 



basis. It follows from ( 12 ) that 

2det{Mhpi±f)>^\\e\\lj>^i.f 



20n 



and likewise for 7p(±e). If fi{M) < 6q, then W{p) = 
{e, /, and opposites}, and this concludes the proof. 

Assume now that ^{M) > 9q. Assume also that 
{e,Mf) < 0, and denote g := —e — f. We have W{p) = 
{e, /, (7,and opposites}, and we obtain as announced 



2det(M)7p(±.g) = -(e. A//) - fi{M) > 0o. 



□ 



Let Th,T'j^ be the energies associated to the stencils 



Let p e ft-Z^ and let ei , • • • , Cfc be the consecutive ele- 
ments of Vh{p), in trigonometric order. We define for all 
l<i<k, denoting M := M{p), 

~h( \ _ (ct - e^-i, Mei_i) + (e, ~ e»+i, Me,+i) 
^P^^'>- 4detM 

with the periodic conventions Ck+i := ei, eg := e^. We 
also set 7p = on Z^ \ {ei, • • • , e^}. 
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Lemma 13. For any p S hZ"^ , one has on T? 



|7p -7pl <ehi, 

where Eh is given 
min{dctM((7); qeil}. 



and - 7p| < CoOhv'p, (53) 

in Lemma and Cq^ := 



Proof. The coefficients 7p, 7^, 7^, are all equal to zero 
outside of W'{p). This holds by construction of 7p, and 
by Lemma 
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for 7p , 7^. We may therefore forget about 
the presence of ry^ in ([53]). 

First inequality. Lemma |6] states that |7p — 7p | < £h on 
Z^, which concludes the proof. 

Second inequality. If n{M) > 9h, then Vh{p) — V{p). 
Comparing the definition of 7^ with that of 7^ ( 13 1 we 
observe that 7^ = 7^ on Z^, which concludes the proof in 
this case. 

Assume now that n{M) < 9^, and denote by (e, /) 



a M-reduced basis. Looking at (12) and denoting 6 
2det(M), we find that 



\Sjp{±e) ~ \\f\\l,\ 
Likewise \Sjp{±f) — \\ 



\{e,Mf)\^^iiM)<eh. 



Ml 



< 9h- Assuming that 



(e, M f) < 0, we obtain in addition 

5^p{±{e + /)) = /i(A^) < and 7p(±(e - /)) = 0. 
Combining the definition of 7^ with the description of 



the stencil Vh{p) in Lemma 10 we obtain that 

U-e.Mf) \ ( {f-e,Mf) 



25t^{e) = 



(/ + e,A//) 



(/ + e,M/> 



In any case |57p(e) — ||/|||/| < 6h. The expressions and 
estimates of 7p at the points —e,f,—f are obtained simi- 
larly. Likewise, using Lemma \lO\ 



2S^';ie + f) 



,Mf) + {e,Mf) 



ife + feV,,{p), 
otherwise. 



In any case \S^p{e + /)| < 9h- The expressions and 
estimates of 7^ at the points — (e + /),e — /, — (e — /) 
are similar. Comparing the above estimates of 7p, 7^, 
we obtain that S\jp — jp\ < 29h on {e, /, e + f,e — f, 
and opposites} = W'{p). Observing that S — 2det(A/) 
> 2Cq^, we conclude the proof. □ 

In the next and last lemma, we control the contribution 
to the energy of a stencil W'{p), p E hi?, in terms of 
the contributions to Th of W{p) and of the neighboring 
stencils W{p + he), e S W{p). This leads to an estimate of 
J-f^ in terms of J-^ , which concludes the proof of Theorem 

m 

Lemma 14. J-[^{u) < CiFh{u), for any u e LF'{Q.h), with 
Ci := 17. 



Proof. Consider a grid point p £ hi?, and denote M :— 
M(p). Assume first that /i(A/) < 6*0, so that W{p) C 
W'{p). Consider also an arbitrary g e W{p) \ W{p), and 
observe that g = e + f for some Af-reduced basis (e, /). 

We set p' p + e and Af :— M(p'). Applying the 
second point of Lemma [9] we find that (e, /) is also a 
A<r'-reduced basis. Indeed we have as required 

V(Af) < 400 = < ||A/"i||-i < Ai(A//)2, 



and the assumption on r follows from (|43j) and (28). 
Therefore 

/ e W{p'), and h-\p' -p) = e e W{p'). (54) 
We obtain 

|u(p + .g)-u(p)p (55) 
= |w(p + e + /)-u(p)|2 

< 2{\u{p + e + f)-u{p + e)\^ + \u{p + e) - u{p)\^) 
= 2{\u{p' + /) - uip')\' + \u{p + e) - u{p)\'). 
Denote, for all q E hi? , 

Th{u; q) := V \u{q + hg) - u(g)p. 



q) := 



g£W[q) 

Hq + hg)~u{q)[' 

g£W'(q) 



Using (55), we obtain 

^'h {u;p)~ {u;p) < Gh (u; p) (56) 

where Qh{u]p) is given by 

4J-;,(«;p) + 2E,6 w(p)^h{u]p + g), \i p,{M.{p)) < 9q 
if /i(M(p)) > 9o 

When Fh{u;p') appears in Gh{u;p), with p,p' E /iZ^, p 7^ 



p' , we have h~^{p' — p) E W{p'), see (54). For each p' E 
hl"^, there are thus at most #{W{p')) < 6 points p E hl'^\ 
{p'} such that Th{u;p') appears in Qh{u]p). Summing 
( [56| over p EVth, we thus obtain T'f^{u) — Th{u) < (4 + 2 x 
6)J-'h{u) (the constant could easily be improved), which 
concludes the proof. □ 



4 Numerical experiments 

We compare our scheme AD-LBR with a family of other 
schemes: finite difference, finite elements, and two sche- 
mes from the image processing literature. We begin with 
a quantitative comparison for the discretization of the 
restoration equation, in a synthetic case where the ex- 
act solution is analytically available for reference. The 
second test case is a qualitative comparison of Coherence- 
Enhancing Diffusion (CED), on a real image and the qual- 
ity assessment is by visual inspection. Finally we present a 
3D implementation of CED, based on AD-LBR, for proof 
of feasibility. 
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4.1 The different schemes 

^ AD-LBR: the scheme presented in this work. 

^ Finite Differences (FD). The gradient and the diver- 
gence are discretized using finite differences, the ma- 
trix D is discretized at the ceUs. More precisely the 
gradient is discretized by 



and the divergence is defined as follows: 

div(DVM),,j = 9^(D"9:rM + Ti^'^dyU).^^^ 

with 

(D"a:rU)i-Hl/2,i 

= 2 (^i+l/2,J + l/2 +^1+1/2.^-1/2) (c'2;W)i+l/2,i, 
5a;(D"aa;U)ij 

= (D"5:r'«)i+l/2j - (D"9a,M),_i/2j 
(D^^9:rU)i+l/2j + l/2 

= I'i+l/2j + l/2 2 ((^2;U)i+l/2j + (f^2:'«)i+i/2j + i) , 

dy{Ty^'^dxU)i^j 

= \ ((D^^9a;U)i+i/2,j + i/2 - (D^^9^u)j+i/2j-l/2 
+ (D2^9j.u),_1/2j + i/2 - (D^^aa,u)i_i/2j-l/2) , 

and similar terms involve dyU. This amounts to dis- 
cretize the term div(DVu) using a stencil of 9 points, 
where the coefficients depend on the matrix D. 

2^ Bilinear Finite Elements (Ql). Bilinear finite ele- 
ments, also referred to as Ql finite elements, are linear 
with respect to each space direction. This amounts 
to use a 9 points stencil, where the coefficients are 
different from the previous scheme. 

^ Scharr-Weickert scheme (SW). This scheme, intro- 
duced in [22], is based on a second order approxima- 
tion of the gradient using a 3 x 3 centered stencil. As 
a result, it offers good accuracy and rotation invari- 
ance when applied to a sufficiently smooth function, 
but lacks robustness guarantees such as the discrete 
maximum principle and spectral correctness (see SjlJ, 
even for D = Id. The stencil for this scheme has size 
5x5. 

2^ Weickert's Non-Negative scheme (W-NN). The co- 
efficients of this scheme, detailed in ^3] page 95, 
are non-negative as long as the anisotropy satisfies 
K.< l-hV2~2.41. 



To fix the ideas and illustrate the difference between 
the schemes, we propose to compute the stencil and the 
coefficients for different constant diffusion tensors D, in 
isotropic and anisotropic cases. Denoting by R the matrix 
of rotation by the angle 9 = 7r/6, and by k > 1 the chosen 
anisotropy ratio, we set, identically on M^: 



D 



R 



1 





R' 



(57) 



The results are presented in Tables [T] and [2] Note that 
for the two last cases (anisotropy k — \/T0 and k — -\/50) 
the AD-LBR stencil contains points that are further than 
the 3x3 neighborhood of the pixel. However the sten- 
cil contains 6 points, as expected. This contrasts with the 
schemes FD, Ql, W-NN where only the 3x3 neighborhood 
is involved. Another observation is that the off-center 
stencil coefficients of the AD-LBR are non-positive (this 
gives non-positive off-diagonal coefficients for — div(DV)), 
in contrast with schemes FD, Ql, SW, and with scheme 
W-NN for anisotropy k > 1 + \f2. This is an essential 
property of AD-LBR, see ([5]), and as a consequence our 
scheme satisfies, unconditionally, the discrete maximum 
principle [1 , 3j . 

The largest eigenvalue of the discrete operator 
— div(DV) is given in Table [sj for the different schemes. 
It turns out that AD-LBR has in most cases the small- 
est eigenvalues among all schemes, excepted scheme SW. 
This property allows (although this was not done in our 
numerical experiments) to use larger time steps for AD- 
LBR than for the other schemes, when solving parabolic 
equations ([2| or ( 60 ) with an explicit time discretization. 



Remark 2 (Additive Operator Splitting). Denote by A 
the symmetric NxN non-negative matrix associated to the 
AD-LBR energy £h ([3]), where N := #{Vlh). Let also U G 
he the vector of values of u G L^(il/i), so that U"^ AU = 
£h{u). The explicit numerical scheme (used in ^4-4 <^i^d, 
P75| ) for the parabolic FDE ^ has the form U' = {I - 
tA)U , where the time step should satisfy < r < 2/||A|| 
for stability. Additive Operator Splitting (AOS) ^ISl allows 
larger time steps, by introducing a semi implicit scheme 



U' 



- y il + mTA,)-^{l~TA)U, 



where A — Aq + ■ ■ ■ + Am + A. For this to be useful, the 
matrices 1 -\- mrAi should be easy to invert. The matrix 
A should be symmetric, non-negative, and its norm \\A\\ 
should be significantly smaller than \\A\\, so that larger 
time steps < r < 2/||y4|| can be used. 

The AD-LBR is completely compatible with this ap- 
proach, since the discrete energy ^ can be written as the 
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sum of the simpler energies: 

+ |u(z — he) ~ u(z)p) , 



where e ranges over elements of 1? greater than (0, 0) in 
lexicographic order, and with co-prime coordinates. The 
N X N symmetric non-negative matrix A^, representing 
the energy is tridiagonal if flh is enumerated in lines 
directed by e; hence the inversion o/ 1 + mrAe has cost 
0{N). AOS could be applied to the decomposition A — 



A,, 



A,, 



A, with suitably chosen ei 



e.g. 



(1,0), (0,1), (1,1), (1,-1). This acceleration procedure 
was not used in our numerical experiments. 



Table 1: The stencil coefficients for different constant dif- 
fusion tensors, and the different schemes presented. The 
value of the anisotropy k is given in the second row, and 
the orientation of the principal axis is = 7r/6, see (57 1. 
The bold coefficient indicates the center node. 



In some 

examples we present for clarity reasons only half of the 
stencil (the other half can be deduced by symmetry) . Sten- 
cil entries are highlighted when they are positive and off- 
center - an undesirable property which gives rise to stabil- 
ity issues. 



K 


K = 1 (D = Id) 




stencil for 
AD-LBR 


0-10 
-1 4 -1 
0-10 


-0.41 -0.22 
-0.66 2.57 -0.66 
-0.22 -0.41 


stencil for 
FD 


0-10 
-1 4 -1 
0-10 


0.11 -0.63 -0.11 
-0.88 3 -0.88 
-0.11 -0.63 0.11 


stencil for 
Ql 


if-' -1 -M 

- -1 8 -1 




-0.14 -0.13 -0.36 
-0.38 2 -0.38 
-0.36 -0.13 -0.14 


stencil for 
SW 


-0.1 -0.06 -0.02 
0.12 -0.06 
0.46 0.12 -0.1 
0.12 -0.06 
-0.1 -0.06 -0.02 


-0.06 -0.05 -0.02 
0.01 -0.04 -0.06 
0.35 0.07 -0.09 
0.01 0.04 -0.04 

-0.06 -0.02 -0.01 


stencil for 
W-NN 


0-10 
-1 4 -1 
0-10 


-0.41 -0.22 
-0.66 2.57 -0.66 
-0.22 -0.41 



4.2 A test case with an explicit solution 

Consider an image v G L'^{il), defined on a domain il, and 
a diffusion tensor field D : — >■ S*^. A classical approach 
to restore the image v, if it has been corrupted by additive 
noise, is to find u G H^(fl) which minimizes: 



+ A 



iVul 



|2 

Id- 



(58) 



Table 2: The stencil coefficients for different metrics and 
the different schemes presented, similarly to Table [T] but 
with more pronounced anisotropics. 



K 






stencil for 
AD-LBR 


-0.26 -0.06 
1.16 -0.26 


-0.11 -0.16 
0.55 -0.01 


stencil for 
FD 


0.19 -0.32 -0.19 
-0.77 2.2 -0.77 
-0.19 -0.32 0.19 


0.21 -0.27 -0.21 
-0.76 2.04 -0.76 
-0.21 -0.27 0.21 


stencil for 
Ql 


0.01 0.04 -0.38 
-0.41 1.47 -0.41 
-0.38 0.04 0.01 


0.04 0.08 -0.38 
-0.42 1.36 -0.42 
-0.38 0.08 0.04 


stencil for 
SW 


-0.02 -0.04 -0.02 
0.09 -0.08 -0.07 
0.25 0.04 -0.08 
0.09 0.08 -0.02 

-0.02 0.004 -0.003 


-0.02 -0.04 -0.02 
0.09 -0.08 -0.07 
0.24 0.03 -0.08 
0.09 0.08 -0.02 

-0.02 0.01 -0.002 


stencil for 
W-NN 


0.06 -0.39 
-0.39 1.42 -0.39 
-0.39 0.06 


0.16 -0.42 
-0.33 1.19 -0.33 
-0.42 0.16 



Table 3: Largest eigenvalue of the discretized operator 
— div(DV), for the constant metric JD — D, where the 
matrix D is given on Tables [l] an d [5] The time step, in 
the explicit discretization of (60), should not exceed the 



inv erse of this value, see Remark [2| 





K = 1 






K = V50 


eigenvalue 
AD-LBR 


8 


4.27 


2.06 


1.06 


eigenvalue 
FD 


8 


6.22 


5.06 


4.85 


eigenvalue 
Ql 


5.7 


4.94 


4.32 


4.20 


eigenvalue 
SW 


1 


1 


1 


1 


eigenvalue 
W-NN 


8 


4.27 


3.1 


3.02 



In other words, u is a penalized least squares approxima- 
tion of V. The parameter A > should be adjusted so as to 
avoid excessive smoothing (for large A), or insufficient de- 
noising (for small A). The solution u can be characterized 
as the solution to the static elliptic PDE: 



-Adiv(DVM) 
Vu.n = 0, 



on fl. 
on dfl. 



(59) 



In applications [TTJ [T3] the diffusion tensor D is usually 
adapted to the local image structure, in order to avoid 
smoothing the edges of v. We construct below a test case 
(image v and tensor field D), for which the solution u is 
known analytically. 

In order to obtain an analytic solution, we first consider 
a separable problem where the image is invariant by trans- 
lation along the horizontal axis, and the metric is constant 
with axes parallel to the coordinate axes. This first prob- 
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lem is invariant under translations along the s-axis, and 
therefore boils down to a 1-dimensional problem. This sep- 
arable problem is then transported by a diffeomorphism 
in order to obtain a new problem where the axes of the 
metric are no more parallel to the coordinate axes. 

The analytical image is composed of a black and a white 
stripe: vo{x) — lxi>o.5i where x = (xi,X2), see Figure [5] 
Given k > 1, we consider the constant diffusion tensor 



Dn 







The analytical solution uq of (59), applied to Dq, vq, 
is known in the case of the infinite domain = M^. In 
Fourier domain all the coefficients are real and: 

This separable problem is transformed using the following 
diffeomorphism: for x = {xi,X2) G ^ 

f{xi,X2) = {xi,X2 + Q!COs(27ra;i)). 

The jacobian of / is 



J[xi,X2) 



-2a sin(27ra;i 



1 



We apply the different restoration schemes to the image 
V — vq o f ^ and the following diffusion tensor: 



D(x) = \ detJix)\.Jix)-^T>oJ{x) 



Jixy^BoJix) 



1/^2 



where we denoted s — 27ra sin(27ra;i). The numerical so- 
lution is compared to the analytical function u = uq o f, 
which is the exact solution in the case of the infinite do- 
main ri = M^. This numerical solution was obtained on 
the bounded domain £7 = [0, Ip, equipped with reflecting 
boundary conditions. Numerical evidence suggests that 
this change of domain and of boundary conditions has 



only an anecdotic impact on the solution of ( 59 1 , with 
the parameters chosen in this test case. 

We used a := 1/3 in the numerical experiments. The 
maximum value of K(D(a;)), among all a; G fi, is equivalent 
to Kmax '■= \/ 1 + {2i:a)'^K ~ 2.3k. 

4.3 Results for the synthetic test case 

We present in Figure |6] the performance results of the dif- 
ferent schemes, for different values of the anisotropy k, ob- 
tained on a series of grids of size ranging from 50 x 50 to 
1200 X 1200. The anisotropy varies from k = 2 to k = 10, 
which are relevant values for imaging applications, see the 
numerical experiments in \ 4.4 The quality of a scheme 





Figure 5: Left: image fana- Right: image v 
transformed by the diffeomorphism /. 



is measured by the difference and the semi-norm 
difference between the numerical solution and the analyt- 
ical solution. Note that the error is concentrated close 
to the discontinuity, since the solution tends rapidly to a 
constant (0 or 1) far from the discontinuity. We chose the 
smoothing parameter A — 10^'^ in (58 1. The linear equa- 



tion obtained by the discretization of ( 59 ) is solved using 
Conjugate Gradient. 

We also tested extreme anisotropics, k > 100 (thus 
'^max > 230), which can be relevant in physics related 
applications. None of the tested schemes showed convinc- 
ing results, and we thus refer to [4] for a radically different 
approach tailored for this setting. 

The performance advantage of the AD-LBR is particu- 
larly clear when the error is measured in the semi- 
norm: for the anisotropy k = 10 and the resolution 
500 X 500, which are relevant values in image processing, 
AD-LBR outperforms its alternatives by a factor ranging 
from 3 to 5. 

4.4 Coherence-enhancing diffusion 

In order to document the interest of our discretization, we 
implement Coherence-Enhancing Diffusion |14j using the 
different numerical schemes at our disposal. The following 
parabolic equation is considered: 



dtu = div(D( Jp(Vw„))V'u). 



(60) 



This equation is non-linear since the diffusion tensor de- 
pends on the solution u, in addition to the four user de- 
fined parameters cr, p, C € M+, a e]0,l[. Let K„ (resp. 
Kp), be the Gaussian kernel of variance a (resp. p). De- 
fine the convolution u^, := *w, and the structure ten- 
sor Jp := Kp -k (Vuo-VuJ). The diffusion tensor D(Jp) 
possesses the same eigenvectors {vi,V2) as Jp, and if the 
eigenvalues of Jp are pi > p2 then the eigenvalues of D(Jp) 
are 



Xi = a 

A2 = a + (1 — a) exp 



-C 
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Comparison k=2 ; c^O.01 



Comparison HI semi norm; k-2 ; c^O.01 




Comparison k=5 ; c^O.01 



Comparison HI semi norm; k-5 ; c^O.01 





Comparison k^lO ; c^O.01 



Comparison HI semi norm; k=1 ; c^O.01 





Figure 6: Left column: relative error in norm (log-log 
scale) for different values of the anisotropy factor: k = 
2, 5, 10; right column: relative error in norm (log-log 
scale) for the same test cases. 



This ensures that one smoothes preferably along the co- 
herence direction V2 , with a diffusivity that increases with 
respect to the coherence (/ii— /i2)^. When the time param- 
eter t becomes large, the image tends to a constant image, 
therefore it is necessary to stop the process at some fi- 
nite time T. The ratio of the eigenvalues is bounded by 
<l/a, hence k < 1/y/a. 



We used an explicit time discretization for (60|), with a 
time step A^. The image at time step (n 

defined by the explicit equation: 



l)At is 



At 



div(D(J,(VO)Vu"). 



The parameters used in our simulation were: a — 0.5, 
p = 4, C = 10-^ a = 10^2 and At = 0.02. This gives 
a maximum anisotropy of k = 10. The algorithm was 
applied to a fingerprint image. The results obtained for 
T = 10 are shown in Figures[7]|8]and[9]and they document 
the ability of our scheme to close interrupted lines more ef- 
ficiently than the other schemes. The largest eigenvalue of 
the discrete operator — div(DV) at t — is given in Table 
|4]for the different schemes. As it was already noticed in the 
constant metric case, it turns out that AD-LBR has the 
smallest eigenvalues among all schemes, excepted scheme 
SW. This property allows (although this was not done in 



Table 4: Largest eigenvalue of the discretized operator 
- div(DV), where D = D( Jp(Vu^)) at i = 0. 



scheme 


AD-LBR 


FD 


Ql 


SW 


W-NN 


eigenvalue 


3.75 


5.67 


5.09 


0.96 


3.83 



our numerical experiments) to use larger time steps for 
AD-LBR than for the other schemes. 

Note also that ridges are clearer, and valleys are darker, 
using AD-LBR than with the other schemes. (Gray-scale 
range is the same for all images, see also Figure 10). This 



reflects the fact that AD-LBR avoids, better than the 
other schemes, smoothing transversally to the orientation 



encoded in the continuous anisotropic PDE ( 60 ) 



Remark 3 (Computation time). Numerical solvers of the 
parabolic PDE (60) combine three main components: (i) 



Constructing the diffusion tensor, (ii) Assembling the dis- 
cretization stencils and the operator sparse matrix. (Hi) 
Performing an explicit time step. Components (i) and 
(ii) are executed exactly the same number of times, while 
step (Hi) is generally more frequent: in order to save CPU 
time, one typically does not update the operator at each 
time step. We produced a C++ implementation of AD- 
LBR, within the Insight Toolkit open source library. Al- 
though our code is neither parallel nor aggressively opti- 
mized, we believe that comparing the CPU times for steps 
(i), (ii) and (Hi) is informative, and allows to estimate the 
additional cost of the AD-LBR which is essentially con- 
tained in step (ii). 

For our 2D Coherence- Enhancing Diffusion ( CED ) ex- 
periment, on the 512 x 512 fingerprint image, (i) takes 
0.21s, (ii) 0.029s, (Hi) 0.007s. For our 3D CED Experi- 
ment, on 100 X 100 X 100 synthetic data, (i) takes L35s, 
(H) 0.55s, (Hi) 0.052s. In both cases, the AD-LBR specific 
step (ii) is dominated by the construction of the diffusion 
tensor (i). Step (ii) is also dominated by the mere cost (Hi) 
of iterations, provided the operator is updated less than 
once every 5 explicit steps in 2D (10 in 3D). To our eyes, 
the limited additional cost (ii) of AD-LBR is acceptable in 
view of the strong theoretical guarantees, and qualitative 
improvements, brought by this scheme. 



4.5 A 3-dimensional experiment 

In order to illustrate the feasibility of our scheme in 3D 
space, we propose a synthetic example of plane-enhancing 
diffusion. The original image is a radially varying image in 
the cube [0, l]'^, and the gray-level of the point x is defined 
to be 

u^{x) = cos (2(r/i?)3) , 

where r = \x\ and i? = 1/2 is a characteristic size of 
the oscillations. This image presents a series of concentric 
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Original image 



AD-LBR ;T=10 



zoom original image 






zoom AD-LBR ; T-10 



zoom FD ; T-10 





Figure 7: From top to bottom and from left to right: Orig- 
inal image (with two regions highlighted); diffused image 
using AD-LBR; FD; Ql; SW; W-NN. Here T = 10. 



Figure 8: Detail of the region on the right. From top to 
bottom and from left to right: original image; diffused 
image using AD-LBR; FD; Ql; SW; W-NN. 



level-sets. We present in Figure 11 the level sets — 0}, Coiiclusioil 



and a slice through the plane z = 0.7. 

The problem is discretized using 100'^ voxels. The image 



is perturbed by 



where n is an additive Gaussian noise of variance a = 
0.5. The reconstructed image is obtained using a 3D 
Coherence-Enhancing Diffusion PDE, similar to the 2D 
one in section IT4l 



dtu ^ diY{T>{Jp{Vu,))Vu), 

where Jp is the structure tensor defined by Jp — Kp * 
{VucrVuJ), Uo- = K„ -ku. The tensor D(Jp) possesses the 
same eigenvectors (wi,W2,^'3) as Jp, and if the eigenvalues 
of Jp are > ^2 > Ma then the eigenvalues of D( Jp) are 



Ai 
A2 

A, 



a 



a - 



(1 
(1 



a) exp 



a) exp 



-C 



(Mi 



-M2)2 

-C 



(mi - M3)^ 

We used the values a = 0.5, p = 4. 
the noisy image u (levelset and 
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where a — 10^^ 
We present in Figure 
planar slice) and the result after 20 time-steps of At 
10-3. 



We introduced in this paper a new numerical scheme, AD- 
LBR, for anisotropic diffusion in image processing. This 
scheme is non-negative, and its stencils have a limited sup- 
port: 6 points in 2D, 14 points in 3D. The former property 
implies that our scheme respects the maximum principle of 
Alvarez, Guichard, Lions and Morel, which is an essential 
feature of parabolic PDEs. 

AD-LBR outperformed all alternatives known to us in 
a quantitative numerical experiment: a test case in which 
approximate numerical solutions are compared against a 
known analytical solution. In a second qualitative test 
case, different schemes were used to enhance a fingerprint 
image. Our scheme appears here to close more efhciently 
the lines of the fingerprint, and to diffuse less orthogo- 
nally to the lines. The is precisely the purpose of the im- 
plemented PDE, coherence enhancing diffusion. We also 
presented a 3-dimensional implementation as a proof of 
feasibility. 

The construction of the stencils of the AD-LBR is 
both original and non-trivial. The computational load 
for this aspect of the algorithm is fortunately not dom- 
inant, thanks to the use of a tool from discrete geometry: 
lattice basis reduction. The AD-LBR also allows to use 
larger time steps than most of its counterparts, in explicit 
discretizations of parabolic equations. 
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Figure 9: Detail of the region on the left. From top to 
bottom and from left to right: original image; diffused 

image using AD-LBR; FD; Ql; SW; W-NN Figure 10: Evolution under CED of a section of the finger- 

print image. The ridges in the evolved image are higher, 
and the valleys are deeper, with AD-LBR than with the 
other schemes. This illustrates the fact that AD-LBR, 
respecting the continuous PDF, diffuses more along the 
structure and less in the orthogonal direction. From top 
to bottom: location of the section of the image; section at 
T = 10; section at T = 50; section at T = 100. 
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